---
title: "2 - Meta analysis"
author: "Alvaro Passi-Solar"
date: "`r Sys.Date()`"
output:
  word_document: default
  html_document: default
---

# Libraries

```{r, include=F}
# birtcohortx<-""
# condx<- "PTSD1"
library(tidyverse)
library(RColorBrewer)
library(weightr)
library(meta)
library(metasens)
library(metafor)
library(captioner)
fig_nums<-captioner(prefix = "Figure 1.",  auto_space = FALSE)
tab_nums<-captioner(prefix = "Table 1.",  auto_space = FALSE)

settings.meta(digits = 3)
options(scipen=999)

```

# Read data from 1 - Preparation_base_SR_Rmd

```{r, echo=FALSE, warning=FALSE, include=F}
df0<-read.csv("df0_SR.csv")


```

# Functions

```{r}

myTryCatch <- function(expr) {
  warn <- err <- NULL
  value <- withCallingHandlers(
    tryCatch(expr, error=function(e) {
      err <<- e
      NULL
    }), warning=function(w) {
      warn <<- w
      invokeRestart("muffleWarning")
    })
  list(value=value, warning=warn, error=err)
}        
```

```{r, echo=FALSE, warning=FALSE, include=F}


FitFlextableToPage <- function(ft, pgwidth = 7){
  
  ft_out <- ft %>% flextable::autofit()
  
  ft_out <- flextable::width(ft_out, width = dim(ft_out)$widths*pgwidth /(flextable::flextable_dim(ft_out)$widths))
  return(ft_out)
}
```







# identify list_any

```{r}

list_any<-c(
  "ANY1", "ANY2", "ANY4",
  "PTSD1", "PTSD2","PTSD4"
  
)
```



# Subset obs based on _birtcohort study

```{r}
# birtcohortx<-"_birtcohort"
# birtcohortx<-""

if(birtcohortx=="_birtcohort"){
  df0_exc<-df0%>%
    subset(!is.na(vi))%>%
    subset(vi>0) %>%
    subset(grepl("Prevalence of comorbidities between mood and anxiety disorders",Title) | 
             grepl("Anxiety disorders in young people",Title) | grepl("Mental disorders and suicide risk in emerging",Title)| 
             grepl("Breastfeeding and mental health in adulthood",Title)| 
             grepl("Posttraumatic stress disorder: Prevalences",Title))
  
  df0<-df0%>%
    subset(!is.na(vi))%>%
    subset(vi>0) %>%
    subset(!grepl("Prevalence of comorbidities between mood and anxiety disorders",Title)) %>%
    subset(!grepl("Anxiety disorders in young people",Title)) %>%
    subset(!grepl("Mental disorders and suicide risk in emerging",Title)) %>%
    subset(!grepl("Breastfeeding and mental health in adulthood",Title)) %>%
    subset(!grepl("Posttraumatic stress disorder: Prevalences",Title))
}
```



```{r, include=F}
# condx<-"AG1"
# condx<-"ANY2"

condx1<-substr(condx,nchar(condx),nchar(condx))

condx0<-substr(condx,1,nchar(condx)-1)

if(condx1=="1"){
  df0_ANY1<-df0%>%
    subset(MAALL_j==paste0(condx0,"1"))
}

if(condx1=="2"){
  df0_ANY1<-df0%>%
    subset(MAALL_j==paste0(condx0,"2"))
}

if(condx1=="4"){
  df0_ANY1<-df0%>%
    subset(MAALL_j==paste0(condx0,"4"))
}

df0_ANY1$ID_original<-df0_ANY1$ID

```




# `r condx`

```{r , include=F}
# "condx" is defined in a loop in 0-Render.Rmd

```

\newpage

# Bias assessment


```{r , include=F, echo=FALSE,  warning = FALSE, out.width="100%", fig.height=20, fig.width=15, fig.align='center'}
m3=rma(yi,vi,
       data=df0_ANY1,
       method="REML",
       weighted=TRUE,
       slab=ID,
       random = ~ 1 | ID,
       digits=5)


summary(m3)
preds <- predict(m3, transf=transf.ilogit)
preds
```


**`r paste0(fig_nums(name = "fpo", caption =""),"Forest plot original")`**

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center'}
forest(m3, 
       slab=paste0(Country,": ",ID_AP),
       atransf=transf.ilogit, 
       refline = predict(m3)$pred,
       order=Country,
       digits=5)


```





```{r, include=F}

m4=rma(yi,vi,
       data=df0_ANY1,
       method="REML",
       weighted=TRUE,
       slab=ID,
       random = ~ 1 | ID,
       digits=5)

pes=predict(m4,
            transf=transf.ilogit
)

# 
summary(m4)
confint(m4)

print(m4,digits=6)
```

\newpage

**`r paste0(fig_nums(name = "doip", caption =""),"lfk-index & doiplot")`**



```{r, include=T}
# The Doi plot substitutes the traditional funnel plot, which displays precision against effect, with a folded normal quantile (Z-score) versus effect plot. When paired with the Galbraith plot, the Doi-Galbraith plot provides a unified graphical summary, offering an immediate visual assessment of both study symmetry and heterogeneity.

lfk.m <- lfkindex(yi, sqrt(vi), data = df0_ANY1)
lfk.m$lfkindex
doiplot(lfk.m,  xlab = "Logit transformed proportion")



```


```{r, include=T}
if (condx %in% list_any & birtcohortx==""){
  # print(condx)
  pdf(file=paste0("Render/",condx,"_doip_",birtcohortx,".pdf")) 
  
  doiplot(lfk.m,lfkindex = FALSE, xlab = "Logit transformed proportion")
  dev.off() 
}


```

## Funnel plot"


\newpage

**`r paste0(fig_nums(name = "funnel", caption =""),"Funnel plot")`**

```{r, include=T}

funnel(m4, level=c(95),atransf=transf.ilogit, legend=FALSE,
       col = "violet",
       back="white", shade="white", hlines="white")



if (condx %in% list_any & birtcohortx==""){
  # print(condx)
  pdf(file=paste0("Render/",condx,"_funnel_",birtcohortx,".pdf")) 
  
  funnel(m4, level=c(95),atransf=transf.ilogit, legend=FALSE,
         col = "violet",
         back="white", shade="white", hlines="white")
  
  dev.off() # Turn the PDF device off
}
```




```{r}
regtest(m4,model="lm", predictor="sei")

```



```{r, include=T, warning=F}

regtest(m4,model="rma", predictor="sei")
```


```{r, include=T, warning=F}

m4$data$vi_inv<-1/m4$data$vi
m4$data$SE_inv<-1/m4$data$SE


m4$data %>%
  ggplot(aes(y=yi, x=vi_inv)) +
  geom_point()+
  geom_smooth(method = "lm", se = FALSE)

ranktest(m4)

fsn(yi,vi,data=df0_ANY1)
```
## Funnel plot trimfill

\newpage

**`r paste0(fig_nums(name = "trimfill", caption =""),"Funnel plot trimfill")`**

```{r, include=T, warning=F}
taf <- trimfill(m4)
# Original
predict1<-predict(m4,
                  transf=transf.ilogit
                  
)
predict1d<-as.data.frame(predict1)
predict1d$model<-"Original"
# Trimfill
predict2<-predict(taf,
                  transf=transf.ilogit)
predict2d<-as.data.frame(predict2)
predict2d$model<-"Trim-fill"
predict12d<-bind_rows(predict1d,predict2d)

```

\newpage

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}

ft1<-flextable::flextable(predict12d)
FitFlextableToPage(ft1)

```


```{r, include=T, warning=F}
### draw funnel plot with missing studies filled in
funnel(taf,level=c(95), atransf=transf.ilogit,
       col = "violet",
       back="white", shade="white", hlines="white")

if (condx %in% list_any & birtcohortx==""){
  
  pdf(file=paste0("Render/",condx,"_fp_trim_",birtcohortx,".pdf")) 
  funnel(taf,level=c(95), atransf=transf.ilogit,
         col = "violet",
         back="white", shade="white", hlines="white")
  dev.off() 
}

```








```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center'}


pes.da=rma(yi,vi,
           data=df0_ANY1,
           method="REML",
           weighted=TRUE,
           slab=ID,
           digits=5)
pes<-predict(pes.da,transf=transf.ilogit)


stud.res=rstudent(pes.da)
abs.z=abs(stud.res$z)


l1o=leave1out(pes.da)
# l1o

yi=l1o$estimate
vi=l1o$se^2

```

\newpage

**`r paste0(fig_nums(name = "inf", caption =""),"Influential")`**


```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center'}
forest(yi,vi,
       slab=paste0(df0_ANY1$Country,": ",df0_ANY1$ID_AP),
       atransf=transf.ilogit, 
       refline = predict(m3)$pred,
       order=df0_ANY1$Country,
       xlab="Summary proportions leaving out each study",
       digits=5)



```




```{r, include=F}

baujat(pes.da)

df0_ANY1$ID[15]
df0_ANY1$ID[1]
df0_ANY1$ID[18]

inf=influence(pes.da)

inf_Df<-inf$inf%>%
  as.data.frame()%>%
  arrange(-abs(rstudent))
inf_Df$ID_o<-rownames(inf_Df)
inf_Df$ID<-rownames(inf_Df)
# print(inf_Df);
plot(inf)
# inf_Df$ID<-gsub(".1","",inf_Df$ID,fixed=T)
# inf_Df$ID<-gsub(".2","",inf_Df$ID,fixed=T)
inf_Df$ID[1]
df0_ANY1$uno<-1
df0_ANY1<-df0_ANY1 %>%
  group_by(ID) %>%
  mutate(ticker=cumsum(uno),
         n_sum=n()) %>%
  select(ID, ticker,n_sum, everything())

df0_ANY1$ID_inf<-paste0(df0_ANY1$ID,".",df0_ANY1$ticker)
df0_ANY1$ID_inf[df0_ANY1$n_sum==1]<-df0_ANY1$ID[df0_ANY1$n_sum==1]

df0_ANY1<-df0_ANY1 %>%
  select(ID, ticker,n_sum, ID_inf,everything())
```






```{r , include=F, echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center'}


pes.da.noutlier=rma(yi,vi,
                    data=subset(df0_ANY1,ID_inf!=inf_Df$ID[1]),
                    method="REML",
                    weighted=TRUE,
                    slab=ID)
pes.noutlier=predict(pes.da.noutlier,
                     transf=transf.ilogit
                     # ,
                     # targ=list(ni=df0_ANY1$total)
)
print(pes,digits=5)
print(pes.noutlier,digits=5)
```


\newpage

**`r paste0(fig_nums(name = "Diagnostic tests", caption =""),"Forest plot without ",inf_Df$ID[1])`**



```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center'}
forest(pes.da.noutlier, 
       slab=paste0(
         subset(df0_ANY1,ID_inf!=inf_Df$ID[1])$Country,
         ": ",
         subset(df0_ANY1,ID_inf!=inf_Df$ID[1])$ID_AP),
       order=subset(df0_ANY1,ID_inf!=inf_Df$ID[1])$Country,
       transf=transf.ilogit,
       digits=5)


```





\newpage


# Meta analysis by country

```{r}
print("Meta analysis by country")
```


```{r , include=F, echo=T,  warning = FALSE}
pred_m1_f<-function(){
  m1<-rma(yi=yi,vi=vi,
          data=df0_ANY1,
          method="REML",
          slab=ID,
          # mods = ~ subgroup,
          random = ~ 1 | ID,
          digits=5)
  
  
  pred_m1<-predict.rma(m1, transf=transf.ilogit)
  
  pred_m1<-pred_m1 %>%
    as.data.frame()
  
  
  
  
  pred_m1$estimate<-m1$b[1]
  pred_m1$stderror<-m1$se
  pred_m1$tau2s<-round(m1$tau2,2)
  # p$tau2<-res$tau2
  pred_m1$Q<-round(m1$QE,2)
  pred_m1$df<-m1$k - m1$p
  pred_m1$p<-metafor:::.pval(m1$QEp,showeq=TRUE)
  pred_m1$I2<-round(m1$I2,1)
  pred_m1$text<-paste0("Heterogeneity: ",
                       " I²=", pred_m1$I2,"%",
                       ", T²=",pred_m1$tau2s,
                       ", Q=",pred_m1$Q,
                       ", df=",pred_m1$df,
                       " (",
                       "p",pred_m1$p,
                       ")")
  list(m1,pred_m1)
}
pred_m1l<-pred_m1_f()
```




```{r ,  warning = FALSE, include=F}

# subgroupx<-"Country"
# subgroupx<-"Sex"
# subx<-unique(df0_ANY1$subgroup)[1]
# subx<-"Brazil"

sg_f<-function(subgroupx) {
  df0_ANY1$subgroup<-df0_ANY1[[subgroupx]]
  df0_ANY1s<-df0_ANY1 %>%
    subset(!is.na(subgroup))
  
  
  # rma by category
  rma_f<-function(subx){
    print(subx)
    df0_ANY1s_s<-subset(df0_ANY1s,subgroup==subx)
    
    d<-rma(yi=yi,vi=vi,
           data=df0_ANY1s_s,
           digits=5,
    )
    
    p<-predict(d,transf=transf.ilogit
    )
    p$subgroup<-subx
    p<-as.data.frame(p)
    p$estimate<-d$b[1]
    p$stderror<-d$se
    p$tau2s<-round(d$tau2,2)
    p$Q<-round(d$QE,2)
    p$df<-d$k - d$p
    p$p<-metafor:::.pval(d$QEp,showeq=TRUE)
    p$I2<-round(d$I2,1)
    p
  }
  
  
  if(length(unique(df0_ANY1s$subgroup))!=1){
    
    pred_m1<-pred_m1_f()[[2]]
    pred_m1<-pred_m1%>%
      as.data.frame()
    pred_m1$subgroup<-"Total"
    pred_m1$ID_AP<-"RE Total"
    total_tex_pi<-paste0(
      "PI: ",
      format(round(pred_m1$pi.lb*100, digits=2), nsmall = 2),
      "-",
      format(round(pred_m1$pi.ub*100, digits=2), nsmall = 2),"")
    
    pred_m1b<-pred_m1
    pred_m1b$subgroup<-"Total"
    pred_m1b$ID_AP<-""
    pred_m1b$pred<-NA
    pred_m1b$ci.lb<-NA
    pred_m1b$ci.ub<-NA
    
    pred_m1$text<-""
    pred_m1<-bind_rows(pred_m1,pred_m1b)
    
    adf_fl<-lapply(unique(df0_ANY1s$subgroup),
                   rma_f)
    
    adf<-do.call(bind_rows,adf_fl)
    adf<-adf%>%
      select(subgroup,everything())
    # as.character(adf$text)
    adf$text<-paste0("Heterogeneity: ",
                     " I²=", adf$I2,"%",
                     ", T²=",adf$tau2s,
                     ", Q=",adf$Q, #X²
                     # \n
                     ", df=",adf$df,
                     " (",
                     "p",adf$p,
                     ")")
    
    adf$tex_pi<-paste0(
      # "",
      #                format(round(adf$pred*100, digits=2), nsmall = 2),
      "PI: ",
      format(round(adf$pi.lb*100, digits=2), nsmall = 2),
      "-",
      format(round(adf$pi.ub*100, digits=2), nsmall = 2),"")
    
    
    
    res_trycatch<- myTryCatch(
      rma.mv(yi, vi, 
             method="REML",
             mods = ~ subgroup, 
             random = ~ subgroup | ID, 
             slab=ID,
             # struct="ID", 
             data=df0_ANY1s, 
             digits=5)
    )
    if (!is.null(res_trycatch$error)){
      print(res_trycatch$error)
      print("optimizer=Nelder-Mead")
      
      res<-  rma.mv(yi, vi, 
                    method="REML",
                    mods = ~ subgroup, 
                    random = ~ subgroup | ID, 
                    slab=ID,
                    # struct="ID", 
                    data=df0_ANY1s, 
                    digits=5, control=list(optimizer="Nelder-Mead"))
    } else {
      res<-res_trycatch$value
    }
    
    
    
    adf$overall_QM<-res$QM
    adf$overall_df<-res$p-1
    adf$overall_QMp<-res$QMp
    
    
    
    
    
    
    #prepare data forest plot:
    pd<-df0_ANY1s %>%
      tibble()%>%
      arrange(Country, Fieldwordmidpoint)
    
    adf0<-adf %>%
      select(subgroup,text,tex_pi,overall_QM,overall_df,overall_QMp)
    adf0$ID_AP<-""
    
    dd<-adf %>%
      select(subgroup,pred,ci.lb,ci.ub)
    dd$ID_AP<-"RE subgroup"
    dd_adf0<-bind_rows(dd,adf0)
    pd<-bind_rows(pd,dd_adf0)
    pd<-bind_rows(pd,pred_m1)
    pd<-pd %>%
      arrange(ID_AP)
    pd$label_ct<-paste0(
      round(pd$cases,0),
      " (",
      round(pd$total,0),
      ")")
    
    pd$label_ct[pd$ID_AP==""]<-""
    pd$label_ct[pd$ID_AP=="Total"]<-""
    pd$label_ct[pd$ID_AP=="RE subgroup"]<-""
    pd$label_ct[pd$ID_AP=="RE Total"]<-""
    
    pd$label<-paste0(format(round(pd$P*100, digits=2), nsmall = 2),
                     " (",
                     format(round(pd$inf_p*100, digits=2), nsmall = 2),
                     "-",
                     format(round(pd$sup_p*100, digits=2), nsmall = 2),
                     ")")
    
    pd$label[is.na(pd$P)]<-
      paste0(format(round(pd$pred[is.na(pd$P)]*100, digits=2), nsmall = 2),
             " (",
             format(round(pd$ci.lb[is.na(pd$P)]*100, digits=2), nsmall = 2),
             "-",
             format(round(pd$ci.ub[is.na(pd$P)]*100, digits=2), nsmall = 2),
             ")")
    
    pd$label[pd$ID_AP==""]<-""
    
    
    
    pd$size<-1/ sqrt(pd$vi)
    pd$size <- (pd$size/20) / max(pd$size, na.rm = T)
    pd$ID_AP[!is.na(pd$text)]<-pd$text[!is.na(pd$text)]
    
    pd$text_caption[pd$ID_AP!=""]<-unique(pd$text_caption[pd$ID_AP==""])
    pd$P[!is.na(pd$pred)]<-pd$pred[!is.na(pd$pred)]
    pd$inf_p[!is.na(pd$ci.ub)]<-pd$ci.ub[!is.na(pd$ci.ub)]
    pd$sup_p[!is.na(pd$ci.lb)]<-pd$ci.ub[!is.na(pd$ci.lb)]
    pd$pred_t_l<-pred_m1$ci.lb[1]
    pd$pred_t<-pred_m1$pred[1]
    pd$pred_t_u<-pred_m1$ci.ub[1]
    pd$nchar<-nchar(pd$label_ct)
    maxlabelct<-max(pd$nchar)
    pd$nchar_max<-maxlabelct-pd$nchar
    
    
    for(i in 1:length(pd$label_ct)){
      pd$label_ct_t[i]<-paste0(pd$label_ct[i],paste0(rep(" ", times=1+pd$nchar_max[i]), collapse = ""))
    }
    pd$nchar<-nchar(pd$label)
    maxlabelct<-max(pd$nchar)
    pd$nchar_max<-maxlabelct-pd$nchar
    
    pd$label[!is.na(pd$tex_pi)]<-pd$tex_pi[!is.na(pd$tex_pi)]
    
    pd$label[pd$label=="" & pd$subgroup=="Total" & pd$ID_AP!="Total"]<-total_tex_pi
    
    
    for(i in 1:length(pd$label_ct)){
      pd$label_2[i]<-paste0(" ",pd$label[i],paste0(rep(" ", times=pd$nchar_max[i]), collapse = ""))
    }
    # nchar(pd$label_ct_t)
    
    pd$label_concat<-paste0(pd$label_ct_t,pd$label_2)
    pd$label_concat<-gsub("NA","", pd$label_concat)
    pd_sub<-pd %>%
      subset(!is.na(subgroup))%>%
      subset(subgroup!="Total")%>%
      select(subgroup)%>%
      unique()
    pd_sub$order<-0
    pd_sub$Fieldwordmidpoint<-0
    pd_sub$ID_AP<-pd_sub$subgroup
    pd_sub$bold<-1
    
    pd<-bind_rows(pd,pd_sub)
    pd<-pd %>%
      select(order, everything())
    
    
    pd$ID_AP[pd$subgroup=="Total" & pd$ID_AP==""]<-"Total"
    pd$order[pd$ID_AP=="Total"]<-7777
    pd$Fieldwordmidpoint[pd$ID_AP=="Total" ]<-7777
    
    pd$order[grepl("RE s",pd$ID_AP)]<-8888
    pd$order[grepl("Hete",pd$ID_AP)]<-9999
    
    pd$Fieldwordmidpoint[grepl("RE s",pd$ID_AP)]<-8888
    pd$Fieldwordmidpoint[grepl("Hete",pd$ID_AP)]<-9999
    
    pd<-pd %>%
      arrange(subgroup,Fieldwordmidpoint,order)

    pd$order_final<-1:nrow(pd)
    
    pd$nchar<-nchar(pd$ID_AP)
    maxlabelct<-max(pd$nchar)
    pd$nchar_max<-maxlabelct-pd$nchar
    
    for(i in 1:length(pd$ID_AP)){
      pd$ID_AP[i]<-paste(pd$ID_AP[i],paste0(rep(" ",times=pd$nchar_max[i]),collapse = ""))
    }
    pd$label_concat[is.na(pd$label_concat)][1]<-"n (N)        %  (95% CI)"
    pd$label_concat[is.na(pd$label_concat)]<-""
    
    pd$ID_AP<-paste0(pd$ID_AP,pd$label_concat)
    pd$ID_AP<-gsub("( ","(",pd$ID_AP,fixed=T)
    pd$ID_AP<-gsub("- ","-",pd$ID_AP,fixed=T)
    pd$nchar<-nchar(pd$ID_AP)
    maxlabelct<-max(pd$nchar)
    pd$nchar_max<-maxlabelct-pd$nchar
    
    for(i in 1:length(pd$ID_AP)){
      pd$ID_AP[i]<-paste(pd$ID_AP[i],paste0(rep(" ",times=pd$nchar_max[i]),collapse = ""))
    }
    
    
    pd<-pd %>%
      arrange(order_final)%>%
      mutate(ID_AP=factor(ID_AP, levels=unique(ID_AP)) )%>%
      select(order_final,everything())
    
    
    pd<-pd %>%
      select(order_final, everything())%>%
      arrange(order_final)
    pd$y_disc<-factor(pd$order_final)
    pd$y_disc<-fct_rev(pd$y_disc)
    pd<-pd %>%
      select(y_disc, everything())
    
    pd$label_concat<-gsub("NA","",pd$label_concat)
    pd$face<-"plain"
    pd$face[pd$bold==1 & !is.na(pd$bold)]<-"bold"
    pd$face[grepl("Total",pd$ID_AP)]<-"bold"
    
    pd$ID_AP<-fct_rev(pd$ID_AP)
    
    # ))
    pd$text_caption<-paste0("Test of Moderators QM (df=",
                            round(res$p-1,0),
                            ") = ",
                            round(res$QM,4),
                            ", p-val = ",
                            format(round(res$QMp,4), nsmall=4))
    pd<-pd%>%
      select(face, everything())
    
    # windows()
    pd$ID_AP<-sub("         ","",pd$ID_AP,fixed=T)
    
    
    
    p1<-ggplot(pd) +
      geom_vline(aes(xintercept = pred_t_l), col="#3A86FF")+
      geom_vline(aes(xintercept = pred_t), col="#3A86FF")+
      geom_vline(aes(xintercept = pred_t_u), col="#3A86FF")+
      
      
      geom_point(aes(x=P, y=y_disc, size=size), shape=15, col="#FFBE0B")+
      geom_errorbarh(aes(xmax = sup_p, xmin = inf_p, 
                         y=y_disc),  height = .2, col="#FFBE0B")+
      scale_y_discrete(
        breaks= pd$y_disc,
        labels = pd$ID_AP
      )+
      
      geom_point(aes(x=pred, y=y_disc), shape=1, col="#FF006E")+
      geom_errorbarh(aes(xmax = ci.ub, xmin = ci.lb, y=y_disc),  height = .2, col="#FF006E")+
      scale_x_continuous(labels = percent)+
      # coord_cartesian(xlim=c(0,max(pd$sup_p, na.rm = T)*1.4))+
      # theme_f()+
      theme(
        text=element_text(family="mono"),
        axis.ticks.y = element_blank(),
        legend.position ="none",
        # aspect.ratio = 1,
        panel.background = element_rect(fill="white"),
        strip.background = element_blank(),
        strip.placement = "outside",
        strip.text = element_blank(),
        axis.text.y = element_text(size=10, colour = "black", angle = 0, hjust = 0, vjust = 0.5,
                                   face = pd$face),
        axis.text.x = element_text(size=10, colour = "black", angle = 0, hjust = 0.5, vjust = 0.5),
        axis.line.x=element_line(color="black"),
        plot.caption.position = "plot",
        plot.caption = element_text(hjust = 0)
      )+
      labs(
        title=condx,
        caption=unique(pd$text_caption),
        y="",
        x="")
    
    list(p1,pd)
  }
  # table(pd$face)
}
```


```{r ,  warning = FALSE, include=F}
l1<-sg_f("Country")
# subgroupx<-"Country"
```


\newpage

**`r paste0(fig_nums(name = "f1", caption =""),"Forest plot (subgroup country)")`**


```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
p1<-l1[[1]]
p1

ggsave(file=paste0("Render/2 - ",condx,"_forest_plot_",birtcohortx,".pdf"), plot=p1, width=15, height=10, dpi=300)


ggsave(file=paste0("Render/png/2 - ",condx,"_forest_plot_",birtcohortx,".png"), plot=p1, width=15, height=10, dpi=600)

```

\newpage

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
pd1<-l1[[2]]
pd1<-pd1%>%
  select(subgroup,ID,ID_AP,label_ct_t,label_2)
ft1<-flextable::flextable(pd1)
FitFlextableToPage(ft1)

```







# Meta analysis by rurality 


```{r ,  warning = FALSE, include=F}
df0_ANY1$Setting<-gsub("Urbano","Urban", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("urbano","Urban", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("urban","Urban", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("mix","Mix", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("Urban","1. Urban", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("Rural","2. Rural", df0_ANY1$Setting)
df0_ANY1$Setting<-gsub("Mix","3. Mix", df0_ANY1$Setting)
```


```{r ,  warning = FALSE, include=F}
l2<-sg_f("Setting")
# subgroupx<-"Country"
```


```{r ,  warning = FALSE, include=F}
# if(length(unique(df0_ANY1$subgroup))==3){
# subgroupx<-"Setting"
df0_ANY1$subgroup<-df0_ANY1[["Setting"]]
res<-  rma.mv(yi, vi, 
              method="REML",
              mods = ~ factor(subgroup)- 1, 
              random = ~ factor(subgroup) | ID, 
              slab=ID,
              # struct="ID", 
              data=df0_ANY1, 
              digits=5, 
              control=list(optimizer="Nelder-Mead"))

```


```{r ,  warning = FALSE}
summary(res)
```

Multiple comparisions

```{r ,  warning = FALSE, include=F}
# glht(res, linfct=cbind(contrMat(rep(1,3), type="Tukey"))), test=adjusted("none")
t_none<-summary(multcomp::glht(res, linfct=cbind(multcomp::contrMat(rep(1,length(unique(df0_ANY1$subgroup))), type="Tukey"))), test=multcomp::adjusted("none"))


broom::tidy(t_none)
t_none

t_none_m<-broom::tidy(t_none)
t_none_m$adjusment<-"none"

# t_none_m<-t_none$model
t_holm<-summary(multcomp::glht(res, linfct=cbind(multcomp::contrMat(rep(1,length(unique(df0_ANY1$subgroup))), type="Tukey"))), test=multcomp::adjusted("holm"))


broom::tidy(t_holm)
t_holm


t_holm_m<-broom::tidy(t_holm)
t_holm_m$adjusment<-"holm"

names(t_none_m)<-gsub(".","",names(t_none_m), fixed=T)
names(t_holm_m)<-gsub(".","",names(t_holm_m), fixed=T)

t_t_setting<-bind_rows(t_none_m,t_holm_m)
t_t_setting$subgroup<-"Setting"
t_t_setting$condx<-condx

write.csv(t_t_setting,paste0("Render/2 - t_t_setting_i",".csv"),row.names = FALSE)


```



\newpage

**`r paste0(tab_nums(name = "t1_allpaircomp", caption =""),"All Pairwise Comparisons - Settings") `**

1. Urban
2. Rural
3. Mix



```{r}

flextable::flextable(t_t_setting)
```



```{r ,  warning = FALSE, include=F}
t_t_setting0<-read.csv("Render/2 - t_t_setting.csv")
t_t_setting_t<-bind_rows(t_t_setting0,t_t_setting)

t_t_setting_t<-unique(t_t_setting_t)
write.csv(t_t_setting_t,"Render/2 - t_t_setting.csv",row.names = FALSE)

```


\newpage

**`r paste0(fig_nums(name = "f1", caption =""),"Forest plot (subgroup country)")`**


```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
p1<-l2[[1]]
p1

ggsave(file=paste0("Render/2 - ",condx,"_forest_plot_urban_",birtcohortx,".pdf"), plot=p1, width=15, height=10, dpi=300)
ggsave(file=paste0("Render/png/2 - ",condx,"_forest_plot_urban_",birtcohortx,".png"), plot=p1, width=15, height=10, dpi=600)

```

\newpage

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
pd1<-l2[[2]]
pd1<-pd1%>%
  select(subgroup,ID,ID_AP,label_ct_t,label_2)
ft1<-flextable::flextable(pd1)
FitFlextableToPage(ft1)

```

\newpage




# Meta regression "Gini","Homicide_Rate","HDI","GII"



```{r, include=F, warning=F}

m2<-rma.mv(yi, vi, 
           method="ML",
           # mods = ~ Gini, 
           random = ~ 1 | ID, 
           slab=ID,
           struct="DIAG", 
           data=df0_ANY1, 
           digits=5)

fx<-function(varx){
  print(varx)
  df0_ANY1$varx<-df0_ANY1[[varx]]
  # df0_ANY1$Gini
  m3<-rma.mv(yi, vi, 
             method="ML",
             mods = ~ varx, 
             random = ~ 1 | ID, 
             slab=ID,
             struct="DIAG", 
             data=df0_ANY1, 
             digits=5)
  
  
  
  summary(m3)
  
  round(exp(coef(summary(m3))[-1,c("estimate", "ci.lb", "ci.ub")]), 2)
  
  preds <- predict(m3, newmods=c(40:60), transf=transf.ilogit)
  (m2$tau2 - m3$tau2) / m2$tau2
  
  valuesx<-c(min(df0_ANY1$varx, na.rm = T),max(df0_ANY1$varx, na.rm = T))
  preds <- predict(m3, newmods=valuesx, transf=transf.ilogit)
  preds<- preds%>%
    as.data.frame()
  preds$valuesx<-valuesx
  preds$model<-"Predicted"
  df0_ANY1$size<-1/ sqrt(df0_ANY1$vi)
  df0_ANY1$size <- (df0_ANY1$size/10) / max(df0_ANY1$size)
  
  scale_col_sr <- function(...){
    ggplot2:::manual_scale(
      'color', 
      values = setNames(c("#F1BB7B", "#FD6467", "#5B1A18", "#D67236", "#0F0D0E", "#85D4E3", "#8D8680"), unique(df0$Country))
      
    )
  }
  
  
  p<-ggplot() +
    geom_ribbon(aes(x=valuesx, ymin=ci.lb, ymax=ci.ub, fill="95% CI"), data=preds)+
    geom_line(aes(x=valuesx, y=pred, linetype="Linear prediction"), data=preds)+
    geom_point(aes(x=varx, y=P, col=Country, size=size), data=df0_ANY1)+
    scale_col_sr()+
    
    scale_size_continuous(guide = "none")+
    xlab(gsub("_"," ",paste(varx,"")))+
    ylab("Predicted probability")+ 
    
    theme(plot.title = element_text(size=16),
          plot.caption = element_text(size=16),
          panel.background = element_rect(fill='transparent'), #transparent panel bg
          plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
          panel.grid.major = element_blank(), #remove major gridlines
          panel.grid.minor = element_blank(), #remove minor gridlines
          legend.background = element_rect(fill='transparent'), #transparent legend bg
          legend.box.background = element_rect(fill='transparent'),
          # legend.box="vertical", 
          legend.margin=margin(),
          legend.position = "bottom",
          legend.title=element_text(size=16), 
          legend.text=element_text(size=16),
          strip.text.x = element_text(size=16, colour = "black", angle = 0),
          legend.key = element_rect(fill = "white", color = NA),
          axis.ticks = element_blank(),
          axis.text.x = element_text(size=16, colour = "black", angle = 0, hjust = 0.5, vjust = 0.5),
          axis.text.y = element_text(size=16, colour = "black", angle = 0, hjust = 1),
          axis.title = element_text(size=16, colour = "black", angle = 0),
          axis.line = element_line(colour = "grey40", size=1, linetype = "solid"),
          panel.grid.minor.x = element_line(colour="transparent", size=0),
          panel.grid.major.x = element_line(colour="grey95", size=0.1),
          panel.grid.minor.y = element_line(colour="transparent", size=0.1),
          panel.grid.major.y = element_line(colour="transparent", size=0.1),
          # panel.background = element_rect(fill="white"),
          strip.text.y = element_text(size = 12, colour = "black", angle = 0),
          # strip.background = element_rect(colour="white", fill="white"),
          panel.spacing = unit(0.5, "lines"))+
    labs(
      title=condx,
      fill="",
      linetype="",
      col=""
      # x=""
    )+
    guides(color=guide_legend(nrow=2,byrow=TRUE))
  
  mr3<-coef(summary(m3))
  mr3$varx<-varx
  
  list(p,mr3)
}

```

```{r}
df0_ANY1$Death_r_disasters
```

```{r, include=F, warning=F}
listx<-lapply(c("Gini","Homicide_Rate",
                "HDI_n_100","HDI_health_100",
                "HDI_edu_100",
                "HDI_inc_100","GII_100"
                # "Total_n_disasters",
                # "Death_r_disasters"
),fx)
```

```{r}
summary(df0_ANY1$Gini)
summary(df0_ANY1$HDI_n_100)
summary(df0_ANY1$GII_100)

summary(df0_ANY1$Homicide_Rate)

# summary(df0_ANY1$Homicide_Rate)


```


\newpage

**`r paste0(fig_nums(name = "f2", caption =""),"Meta-regression (predicted vs observed)")`**

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
listx[[1]][[1]]
listx[[2]][[1]]
listx[[3]][[1]]
listx[[4]][[1]]
listx[[5]][[1]]
listx[[6]][[1]]
listx[[7]][[1]]
# listx[[8]][[1]]
# listx[[9]][[1]]
list_plots<-c("Gini","Homicide_Rate",
              "HDI_n_100","HDI_health_100",
              "HDI_edu_100",
              "HDI_inc_100","GII_100"
              # "Total_n_disasters",
              # "Death_r_disasters"
)
```


```{r, echo=F, warning=F}


if (condx %in% list_any & birtcohortx==""){
  for(i in 1:length(list_plots)) {
    ggsave(file=paste0("Render/2 - ",condx,"_bubbleplot",list_plots[i],"_",birtcohortx,".pdf"), plot=listx[[i]][[1]], width=15, height=10, dpi=300)
    ggsave(file=paste0("Render/png/2 - ",condx,"_bubbleplot",birtcohortx,".png"), plot=p1, width=15, height=10, dpi=600)
    
  }
}



```


\newpage

**`r paste0(tab_nums(name = "t1", caption =""),"Meta-regression (Odds ratios)")`**


```{r, echo=F, warning=F}
tt<-bind_rows(listx[[1]][[2]],
              listx[[2]][[2]],
              listx[[3]][[2]],
              listx[[4]][[2]],
              listx[[5]][[2]],
              listx[[6]][[2]],
              listx[[7]][[2]]
              # listx[[8]][[2]],
              # listx[[9]][[2]]
)


# tt<-bind_rows(mr1,mr3,mr4,mr5,mr6,mr7,mr8,mr9)
tt$var<-rownames(tt)
rownames(tt)<-NULL
tt<-tt%>%
  subset(!grepl("int",var))%>%
  select(varx, everything())

tt$OR<-round(exp(tt$estimate),3)
tt$OR_l<-round(exp(tt$ci.lb),3)
tt$OR_u<-round(exp(tt$ci.ub),3)
tt$pval<-round(tt$pval,4)


tt<-tt%>%
  subset(!grepl("int",var))%>%
  select(varx, OR,OR_l,OR_u, pval)


flextable::flextable(tt)
```



\newpage

**`r paste0(tab_nums(name = "t12", caption =""),"Meta-regression multivar HDI (coef)")`**


```{r, include=T, echo=F}
m3<-rma.mv(yi, vi, 
           method="ML",
           mods = ~HDI_health_100+ HDI_edu_100+ HDI_inc_100,
           random = ~ 1 | ID, 
           slab=ID,
           # struct="DIAG", 
           data=df0_ANY1, 
           digits=5)
# predict(m2, transf=transf.ilogit)
# summary(m3, transf=transf.ilogit)
bt1<-broom::tidy(m3)
bt2<-round(exp(coef(summary(m3))[-1,c("estimate", "ci.lb", "ci.ub")]), 2)
bt2$var<-rownames(bt2)
bt2<-bt2 %>%
  select(var,everything())
flextable::flextable(bt1)
```


\newpage

**`r paste0(tab_nums(name = "t123", caption =""),"Meta-regression multivar HDI (OR)")`**


```{r, include=T, echo=F}
flextable::flextable(bt2)


```


\newpage

# Meta-analysis by gender

```{r, include=F}
pred_m1_f()

# df0$varx<-as.factor(df0$Sex)
# df0$Outcome_vff



if(condx1=="1"){
  df0_ANY1<-df0%>%
    subset(MASex_j==paste0(condx0,"1")|MASex_j==paste0(condx0,"2"))
}

if(condx1=="2"){
  df0_ANY1<-df0%>%
    subset(MASex_j==paste0(condx0,"3")|MASex_j==paste0(condx0,"4"))
}

if(condx1=="4"){
  df0_ANY1<-df0%>%
    subset(MASex_j==paste0(condx0,"7")|MASex_j==paste0(condx0,"8"))
}





```


```{r ,  warning = FALSE, include=F}
l1<-sg_f("Sex")
```



\newpage

**`r paste0(fig_nums(name = "f1s", caption =""),"Forest plot (subgroup sex)")`**


```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
p2<-l1[[1]]

p2
ggsave(file=paste0("Render/2 - ",condx,"_forest_plot_sex",birtcohortx,".pdf"), plot=p2, width=15, height=10, dpi=300)
ggsave(file=paste0("Render/png/2 - ",condx,"_forest_plot_sex",birtcohortx,".png"), plot=p2, width=15, height=10, dpi=600)
if(is.null(p2)){
  file.remove(paste0("Render/2 - ",condx,"_forest_plot_sex",birtcohortx,".pdf"))
  file.remove(paste0("Render/png/2 - ",condx,"_forest_plot_sex",birtcohortx,".png"))
  
}
```

\newpage


```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
pd2<-l1[[2]]
pd2<-pd2%>%
  select(subgroup,ID,ID_AP,label_ct_t,label_2)
ft1<-flextable::flextable(pd2)
FitFlextableToPage(ft1)

```



```{r, include=F, warning=F}
fx_sex<-function(sexx){
  # sexx<-"F"
  # varx<-"Gini"
  fx<-function(varx){
    print(varx)
    df0_ANY1_sex<-df0_ANY1
    # df0_ANY1_sex$Sex
    df0_ANY1_sex$varx<-df0_ANY1_sex[[varx]]
    df0_ANY1_sex<-df0_ANY1_sex%>%
      subset(Sex==sexx)
    # df0_ANY1$Gini
    m3<-rma.mv(yi, vi, 
               method="ML",
               mods = ~ varx, 
               random = ~ 1 | ID, 
               slab=ID,
               struct="DIAG", 
               data=df0_ANY1_sex, 
               digits=5)
    
    
    
    # predict(m2, transf=transf.ilogit)
    summary(m3)
    
    round(exp(coef(summary(m3))[-1,c("estimate", "ci.lb", "ci.ub")]), 2)
    
    preds <- predict(m3, newmods=c(40:60), transf=transf.ilogit)
    (m2$tau2 - m3$tau2) / m2$tau2
    
    valuesx<-c(min(df0_ANY1_sex$varx, na.rm = T),max(df0_ANY1_sex$varx, na.rm = T))
    preds <- predict(m3, newmods=valuesx, transf=transf.ilogit)
    preds<- preds%>%
      as.data.frame()
    preds$valuesx<-valuesx
    
    df0_ANY1_sex$size<-1/ sqrt(df0_ANY1_sex$vi)
    df0_ANY1_sex$size <- (df0_ANY1_sex$size/10) / max(df0_ANY1_sex$size)
    
    
    p<-ggplot() +
      geom_ribbon(aes(x=valuesx, ymin=ci.lb, ymax=ci.ub, fill="95% CI"), data=preds)+
      scale_fill_manual(values=c("#f5cc46"))+
      geom_line(aes(x=valuesx, y=pred, linetype="Linear prediction"), data=preds)+
      # ggnewscale::new_scale_color() +
      geom_point(aes(x=varx, y=P, col=Country, size=size), data=df0_ANY1_sex)+
      scale_size_continuous(guide = "none")+
      
      xlab(gsub("_"," ",paste(varx,"")))+
      ylab("Predicted probability")+ 
      theme(plot.title = element_text(size=16),
            plot.caption = element_text(size=16),
            panel.background = element_rect(fill='transparent'), #transparent panel bg
            plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
            panel.grid.major = element_blank(), #remove major gridlines
            panel.grid.minor = element_blank(), #remove minor gridlines
            legend.background = element_rect(fill='transparent'), #transparent legend bg
            legend.box.background = element_rect(fill='transparent'),
            # legend.box="vertical", 
            legend.margin=margin(),
            legend.position = "bottom",
            legend.title=element_text(size=16), 
            legend.text=element_text(size=16),
            strip.text.x = element_text(size=16, colour = "black", angle = 0),
            legend.key = element_rect(fill = "white", color = NA),
            axis.ticks = element_blank(),
            axis.text.x = element_text(size=16, colour = "black", angle = 0, hjust = 0.5, vjust = 0.5),
            axis.text.y = element_text(size=16, colour = "black", angle = 0, hjust = 1),
            axis.title = element_text(size=16, colour = "black", angle = 0),
            axis.line = element_line(colour = "grey40", size=1, linetype = "solid"),
            panel.grid.minor.x = element_line(colour="transparent", size=0),
            panel.grid.major.x = element_line(colour="grey95", size=0.1),
            panel.grid.minor.y = element_line(colour="transparent", size=0.1),
            panel.grid.major.y = element_line(colour="transparent", size=0.1),
            # panel.background = element_rect(fill="white"),
            strip.text.y = element_text(size = 12, colour = "black", angle = 0),
            # strip.background = element_rect(colour="white", fill="white"),
            panel.spacing = unit(0.5, "lines"))+
      labs(
        title=condx,
        fill="",
        linetype="",
        col=""
        # x=""
      )+
      guides(color=guide_legend(nrow=2,byrow=TRUE))
    
    mr3<-coef(summary(m3))
    mr3$varx<-varx
    
    list(p,mr3)
  }
  
  lapply(c("Gini","Homicide_Rate",
           "HDI_n_100","HDI_health_100",
           "HDI_edu_100",
           "HDI_inc_100","GII_100"
           # "Total_n_disasters",
           # "Death_r_disasters"
  ),fx)
  
}


listxx<-lapply(c("M","F"),fx_sex)
```


\newpage

**`r paste0(fig_nums(name = "f2_males", caption =""),"Males: Meta-regression (predicted vs observed)")`**

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
listx<-listxx[[1]]
listx[[1]][[1]]
listx[[2]][[1]]
listx[[3]][[1]]
listx[[4]][[1]]
listx[[5]][[1]]
listx[[6]][[1]]
listx[[7]][[1]]
# listx[[8]][[1]]
# listx[[9]][[1]]
```


```{r, echo=F, warning=F}
list_any<-c(
  "ANY1", "ANY2", "ANY4",
  "PTSD1", "PTSD2","PTSD4"
)

if (condx %in% list_any & birtcohortx==""){
  for (i in 1:length(list_plots)) {
    ggsave(file=paste0("Render/2 - ",condx,"_males_bubbleplot",list_plots[i],"_",birtcohortx,".pdf"), plot=listx[[i]][[1]], width=15, height=10, dpi=300)
    ggsave(file=paste0("Render/png/2 - ",condx,"_males_bubbleplot",list_plots[i],"_",birtcohortx,".png"), plot=listx[[i]][[1]], width=15, height=10, dpi=600)
  }
}



```


\newpage

**`r paste0(fig_nums(name = "f2_females", caption =""),"Females: Meta-regression (predicted vs observed)")`**

```{r , echo=FALSE,  warning = FALSE, out.width="100%", fig.height=10, fig.width=15, fig.align='center', include=T}
listx<-listxx[[2]]
listx[[1]][[1]]
listx[[2]][[1]]
listx[[3]][[1]]
listx[[4]][[1]]
listx[[5]][[1]]
listx[[6]][[1]]
listx[[7]][[1]]
# listx[[8]][[1]]
# listx[[9]][[1]]

if (condx %in% list_any & birtcohortx==""){
  for (i in 1:length(list_plots)) {
    ggsave(file=paste0("Render/2 - ",condx,"_females_bubbleplot",list_plots[i],"_",birtcohortx,".pdf"), plot=listx[[i]][[1]], width=15, height=10, dpi=300)
    ggsave(file=paste0("Render/png/2 - ",condx,"_females_bubbleplot",list_plots[i],"_",birtcohortx,".png"), plot=listx[[i]][[1]], width=15, height=10, dpi=600)
  }
}

```

\newpage

**`r paste0(tab_nums(name = "t1_sex", caption =""),"Meta-regression (Odds ratios) by sex")`**


```{r, echo=F, warning=F}
listx<-listxx[[1]]

tt1<-bind_rows(listx[[1]][[2]],
               listx[[2]][[2]],
               listx[[3]][[2]],
               listx[[4]][[2]],
               listx[[5]][[2]],
               listx[[6]][[2]],
               listx[[7]][[2]]
               # listx[[8]][[2]],
               # listx[[9]][[2]]
)
tt1$sex<-"Males"

listx<-listxx[[2]]

tt2<-bind_rows(listx[[1]][[2]],
               listx[[2]][[2]],
               listx[[3]][[2]],
               listx[[4]][[2]],
               listx[[5]][[2]],
               listx[[6]][[2]],
               listx[[7]][[2]]
               # listx[[8]][[2]],
               # listx[[9]][[2]]
)
tt2$sex<-"Females"

tt<-bind_rows(tt1,tt2)

# tt<-bind_rows(mr1,mr3,mr4,mr5,mr6,mr7,mr8,mr9)
tt$var<-rownames(tt)
rownames(tt)<-NULL
tt<-tt%>%
  subset(!grepl("int",var))%>%
  select(sex,varx, everything())

tt$OR<-round(exp(tt$estimate),3)
tt$OR_l<-round(exp(tt$ci.lb),3)
tt$OR_u<-round(exp(tt$ci.ub),3)
tt$pval<-round(tt$pval,4)


tt<-tt%>%
  subset(!grepl("int",var))%>%
  select(sex,varx, OR,OR_l,OR_u, pval)


flextable::flextable(tt)
```



